#ifndef _FUNC_CPP_
#define _FUNC_CPP_

#include<iostream>
#include<fstream>
#include "Polynomial.h"
#include "Polynomial.cpp"
#include "Spline.h"
#include "Spline.cpp"
#include "InterConditions.h"
#include <cmath>

double fe(double x)
{
  return 2.0*(sqrt(x)+sqrt(3-x*x))/3;
}
double fe1(double x)
{
  return 2.0*(sqrt(x)-sqrt(3-x*x))/3;
}

double fc(double x)
{
  return 1.0/(1+x*x);
}

double fb(double x)
{
  return 1.0/(1+25*x*x);
}

void B()
{
  ofstream fout;
  fout.open("P-b.txt");
  int N=6;
  for(;N<85;N++)
    {
      if(N==6 || N==11 || N==21 || N==41 || N==81){
	vector<double> x(N);
	vector<double> y(N+2);
	for(int i=0;i<N;i++)
	  {
	    x[i]=-1+i*2.0/(N-1);
	    y[i]=fb(x[i]);
	  }
	y[N]=50.0/(26*26);y[N+1]=-y[N];
	InterConditions p(N,x,y);
	Spline<1,4,ppForm> S=interpolate<4>(p,complete);
	Spline<1,4,ppForm> S1=interpolate<4>(p,notAknot);
	double maxnorm=0;
	double m;
	for(int i=0;i<N-1;i++)
	  {
	    m=abs(S.s_get_value((x[i]+x[i+1])/2)-fb((x[i]+x[i+1])/2));
	    if(m>maxnorm) maxnorm=m;
	  }
	fout<<"For N= "<<N<<endl;
	fout<<"complete cubic spline,the maxnorm of error vector at the mid-points of subintervals:"<<maxnorm<<endl;
	maxnorm=0;
	for(int i=0;i<N-1;i++)
	  {
	    m=abs(S1.s_get_value((x[i]+x[i+1])/2)-fb((x[i]+x[i+1])/2));
	    if(m>maxnorm) maxnorm=m;
	  }
	fout<<"not-a-knot cubic spline,the maxnorm of error vector at the mid-points of subintervals:"<<maxnorm<<endl;
	fout<<"for complete cubic spline:"<<endl; 
	for(double k=-1;k<=1;k+=0.02)
	  fout<<"( "<<k<<","<<S.s_get_value(k)<<" )"<<endl;
	fout<<"( "<<"1"<<","<<S.s_get_value(1)<<" )"<<endl;
	fout<<"for not-a-knot cubic spline:"<<endl;
	for(double k=-1;k<=1;k+=0.02)
	  fout<<"( "<<k<<","<<S1.s_get_value(k)<<" )"<<endl;
	fout<<"( "<<"1"<<","<<S1.s_get_value(1)<<" )"<<endl;  
      }
    }
  fout<<flush;
  fout.close();
}



void C()
{
  ofstream fout;
  fout.open("P-c.txt");
  vector<double> x1(11);
  vector<double> y1(13);
  for(int i=0;i<11;i++)
    {
      x1[i]=-6+i+1;
      y1[i]=fc(x1[i]);
    }
  y1[11]=10.0/(26*26);y1[12]=-y1[11];
  InterConditions p1(11,x1,y1);
  Spline<1,4,cardinalB> S1=interpolate<4>(p1,cubic);
  fout<<"cubic:4.58:"<<endl;
  double k=-5;
  for(int i=0;k<=5.0;i++)
    {
      fout<<"( "<<k<<","<<S1.s_get_value(k)<<" )"<<endl;
      k+=0.1;
    }
  vector<double> x2(10);
  vector<double> y2(12);
  for(int i=0;i<10;i++)
    {
      x2[i]=i-4.5;
      y2[i]=fc(x2[i]);
    }
  y2[10]=fc(-5.0);y2[11]=fc(5.0);
  InterConditions p2(10,x2,y2);
  Spline<1,3,cardinalB> S2=interpolate<3>(p2,quadratic);
  fout<<"quadratic:4.59:"<<endl;
  double k2=-5;
  for(int i=0;k2<=5.0;i++)
    {
      fout<<"( "<<k2<<","<<S2.s_get_value(k2)<<" )"<<endl;
      k2+=0.1;
    }
  fout<<flush;
  fout.close();
}

void D()
{
  ofstream fout;
  fout.open("P-d.txt");
  vector<double> x1(11);
  vector<double> y1(13);
  for(int i=0;i<11;i++)
    {
      x1[i]=-6+i+1;
      y1[i]=fc(x1[i]);
    }
  y1[11]=10.0/(26*26);y1[12]=-y1[11];
  InterConditions p1(11,x1,y1);
  Spline<1,4,cardinalB> S1=interpolate<4>(p1,cubic);
   vector<double> x2(10);
  vector<double> y2(12);
  for(int i=0;i<10;i++)
    {
      x2[i]=i-4.5;
      y2[i]=fc(x2[i]);
    }
  y2[10]=fc(-5.0);y2[11]=fc(5.0);
  InterConditions p2(10,x2,y2);
  Spline<1,3,cardinalB> S2=interpolate<3>(p2,quadratic);
  double k=-3.5;
  fout<<"Error of cubic:"<<endl;
  for(int i=0;k<=3.5;k+=0.5)
    if(k==-3.5 || k==-3 || k==-0.5||k==0||k==0.5||k==3||k==3.5)
      fout<<"at x="<<k<<": "<<abs(S1.s_get_value(k)-fc(k))<<endl; 
  fout<<"Error of quadratic:"<<endl;
  k=-3.5;
  for(int i=0;k<=3.5;k+=0.5)
    if(k==-3.5 || k==-3 || k==-0.5||k==0||k==0.5||k==3||k==3.5)
      fout<<"at x="<<k<<": "<<abs(S2.s_get_value(k)-fc(k))<<endl;
  fout<<"Some of the errors close to machine precision,because in which the sites are the interpolation sites of corresponding spline,that is,at these sites,the interpolation values are very close to its real values,almost the same.According to the nature of machine operation,at these sites,the $E_{S}(x)$ is close to machine precision."<<endl; 
  fout<<"In addition,obviously,cubic cardinal B-spline is more accurate!"<<endl;
  fout<<flush;
  fout.close();
}
void E()
{
  ofstream fout;
  fout.open("P-e.txt");
  fout<<"For n=10:"<<endl;
  double x[6]={0,0.4,0.83,1.73,0.4,0};
  double y[6]={fe(0),fe(0.4),fe(0.83),fe(1.73),fe1(0.4),fe1(0)};
  vector<Vec<double,2> > p(6);
  for(int i=0;i<6;i++)
    {
      p[i][0]=x[i];
      p[i][1]=y[i];
    }
  Spline<2,4,ppForm> S=fitCurve<4>(p,notAknot);
  double k=0;
  double l=S.get_T();
  for(int i=0;k<=l;)
    {
      fout<<S.s1_get_value(k)<<endl;
      k+=0.05;
    }
  fout<<S.s1_get_value(l)<<endl;
  k-=0.05;
  for(int i=0;k>=0;k-=0.05)
    fout<<"( "<<-S.s1_get_value(k)[0]<<","<<S.s1_get_value(k)[1]<<" )"<<endl;
  fout<<S.s1_get_value(0)<<endl;
  fout<<"*********************************"<<endl;
  fout<<"For n=40:"<<endl;
  double x1[21];
  double y1[21];
  double a=0;
  for(int i=0;i<11;i++)
    {
      x1[i]=a;
      a+=0.173;
      y1[i]=fe(x1[i]);
    }
  for(int i=11;i<21;i++)
    {
      x1[i]=x1[20-i];
      y1[i]=fe1(x1[i]);
    }
  vector<Vec<double,2> > p1(21);
  for(int i=0;i<21;i++)
    {
      p1[i][0]=x1[i];
      p1[i][1]=y1[i];
    }
  Spline<2,4,ppForm> S1=fitCurve<4>(p1,notAknot);
  double k1=0;
  double l1=S1.get_T();
  for(int i=0;k1<=l1;)
    {
      fout<<S1.s1_get_value(k1)<<endl;
      k1+=0.05;
    }
  fout<<S1.s1_get_value(l1)<<endl;
  k1-=0.05;
  for(int i=0;k1>=0;k1-=0.05)
    fout<<"( "<<-S1.s1_get_value(k1)[0]<<","<<S1.s1_get_value(k1)[1]<<" )"<<endl;
  fout<<S1.s1_get_value(0)<<endl;
  fout<<"*********************************"<<endl;
  fout<<"For n=160:"<<endl;
  double x2[81];
  double y2[81];
  double a1=0;
  for(int i=0;i<41;i++)
    {
      x2[i]=a1;
      a1+=0.043;
      y2[i]=fe(x2[i]);
    }
  for(int i=41;i<81;i++)
    {
      x2[i]=x2[80-i];
      y2[i]=fe1(x2[i]);
    }
  vector<Vec<double,2> > p2(81);
  for(int i=0;i<81;i++)
    {
      p2[i][0]=x2[i];
      p2[i][1]=y2[i];
    }
  Spline<2,4,ppForm> S2=fitCurve<4>(p2,notAknot);
  double k2=0;
  double l2=S2.get_T();
  for(int i=0;k2<=l2;)
    {
      fout<<S2.s1_get_value(k2)<<endl;
      k2+=0.05;
    }
  fout<<S2.s1_get_value(l2)<<endl;
  k2-=0.05;
  for(int i=0;k2>=0;k2-=0.05)
    fout<<"( "<<-S2.s1_get_value(k2)[0]<<","<<S2.s1_get_value(k2)[1]<<" )"<<endl;
  fout<<S2.s1_get_value(0)<<endl;
  fout<<flush;
  fout.close();
}


void F()
{
  ofstream fout;
  fout.open("P-f.txt");
  vector<double> x1={1,2,5,6,7,8,10,13,17};
  vector<double> y1={3.0,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5,1,-0.67};
  vector<double> x2={17,20,23,24,25,27,27.7};
  vector<double> y2={4.5,7.0,6.1,5.6,5.8,5.2,4.1,3.0,-4.0};
  vector<double> x3={27.7,28,29,30};
  vector<double> y3={4.1,4.3,4.1,3.0,0.33,-1.5};
  InterConditions p1(9,x1,y1);
  InterConditions p2(7,x2,y2);
  InterConditions p3(4,x3,y3);
  Spline<1,4,ppForm> S1=interpolate<4>(p1,complete);
  Spline<1,4,ppForm> S2=interpolate<4>(p2,complete);
  Spline<1,4,ppForm> S3=interpolate<4>(p3,complete);
  double k=1;
  for(int i=0;k<=17;k+=0.1)
    fout<<"( "<<k<<","<<S1.s_get_value(k)<<" )"<<endl;
  for(int i=0;k<=27.7;k+=0.1)
    fout<<"( "<<k<<","<<S2.s_get_value(k)<<" )"<<endl;
  for(int i=0;k<=30;k+=0.1)
    fout<<"( "<<k<<","<<S3.s_get_value(k)<<" )"<<endl;
  fout<<flush;
  fout.close();
}

// void G()
// {
  
// }

#endif
